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We describe a new statistical pattern in the chaotic dynamics of closed inflationary cosmologies, 
associated with the partition of the Hamiltonian into rotational motion energy and hyperbolic 
motion energy pieces, in a linear neighborhood of the saddle-center present in the phase space of the 
models. The hyperbolic energy of orbits visiting a neighborhood of the saddle-center has a random 
distribution with respect to the ensemble of initial conditions, but the associated histograms define 
a statistical distribution law of the form p{x) — C x~'' , for almost the whole range of hyperbolic 
energies considered. We present numerical evidence that 7 determines the dimension of the fractal 
basin boundaries in the ensemble of initial conditions. This distribution is universal in the sense 
that it does not depend on the parameters of the models and is scale invariant. We discuss possible 
physical consequences of this universality for the physics of inflation. 
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The general dynamics of closed Friedmann-Robertson- 
Walker (FRW) inflationary cosmologies, -with a mini- 
mally coupled massive scalar field and perfect fluid or 
-with anisotropy, is chaotic ||,|^. The chaotic behaviour 
has a definite homoclinic origin, due to infinite transver- 
sal crossings of homoclinic cylinders present in the phase 
space of the models, analogous to the case of breaking 
and crossing of homoclinic connections in Poincare's ho- 
moclinic phenomena |^-|^. The topology of homoclinic 
cylinders is a consequence of a critical point of the saddle- 
center type in phase space, and the breaking and 
crossings of the cylinders are due to the non-integrability 
of the dynamics. These homoclinic phenomena produce 
chaotic sets in phase space, in the sense that the dynam- 
ics is highly sensitive to fluctuations in initial conditions 
taken on these sets. The aim of this communication is to 
describe a new statistical pattern in the chaotic dynamics 
engendered by the saddle center and -what are its possible 
consequences in the physics of inflation. Throughout the 
paper we use units such that 8 tt G = 1 = c. 

Our discussion will be basically restricted to closed 
Friedmann-Robertson- Walker (FRW) inflationary cos- 
mologies, with a minimally coupled massive scalar field 
and perfect fluid, but the same construction will apply to 
a Bianchi IX anisotropic cosmology and yields analogous 
results, as we will see. The equations of motion may be 
derived from the Hamiltonian constraint 
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-V{^)-Eo = Q, (1) 



where a and Lp are the scale factor and the scalar field, 
respectively, and pa and are the momenta canonically 
conjugate to a and ip. We consider here a perfect fluid de- 
scribing radiation, and the potential V{ip) = A 



The constant A is considered to correspond to the vac- 
uum energy of the scalar field Lp, and plays the role 
of a cosmological constant term. This model may be 
thought to describe a classical universe just exiting a 
Planck era. Eq is a constant arising from the first in- 
tegral of the Bianchi identities, and proportional to the 
total matter energy of the models. The dynamical sys- 
tem generated by (1) has a saddle-center critical point 
in the finite region of the phase space with coordinates 

E : ip = {),a = aa ^ \f^,Pa = 0,?^ = 0. The energy 
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associated to the critical point E is Eq — Ecrit = 
The phase space also has a pair of critical points at in- 
finity, corresponding to the DeSitter solution, one acting 
as an attractor, and the other as a repeller; the DeSitter 
attractor characterizes an exit to infiation for the models. 

In accord to a theorem by Moser it is always pos- 
sible to find a set of canonical variables such that, in a 
linear neighborhood of a saddle-center, the Hamiltonian 
dynamics is separable into rotational motion and hyper- 
bolic motion pieces, with the associated partial energies 
Erot and Ehyp approximately conserved. The canonical 
variables {a,ip,Pa,Pip) behave as Moser's variables in a 
small neighborhood of E; in fact, expanding the Hamil- 
tonian (1) about E we obtain 

H = Ehyp - Erot + {Ecr - Eo) + 0(3) = 0, (2) 



■a^Sip'^ and Ehyp = J2PI 
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where Erot = + ^' 
6Sa^, with Sa = a — oq, 5ip = (p, and {Ecr — Eq) small. 
0(3) denotes higher order terms which are neglected in a 
small neighborhood of E. In this approximation the scale 
factor a{t) has pure hyperbolic motion, and is completely 
decoupled from the scalar field pure rotational motion in 



1 



the plane {ip,p^p). 

We will now focus our attention to chaotic aspects in 
the dynamics, leading us into a statistical treatment of 
possible experiments made in the models. Our descrip- 
tion will be restricted to sets of orbits which visit a linear 
neighborhood of the saddle-center E, and makes exten- 
sive use of the partition of energy Ehyp — Erot + [Ecrit — 
Eq) w that occurs in this neighborhood. Let us con- 
sider a set V of arbitrary initial conditions which corre- 
spond to initially expanding universes with energy Eq, 
and such that all the orbits visit a linear neighborhood 
of E. The set V may be represented as a small 4-dim vol- 
ume, contained in a 4-dim sphere of radius R (of order 
10~^, 10~^, depending on the order of fluctuations 
we admit in the initial conditions) constructed about the 
point, say, (0.4,0,3.5638181771,0) in phase space, and 
satisfying the Hamiltonian constraint (1). The above 
point belongs to the invariant plane p^p = (p = 0, and 
is located on one of the separatrices reaching the saddle- 
center E. We denote R as the characteristic radius of 
the four volume V. One of the chaotic aspects we sin- 
gle out for this set is the definition of a chaotic exit to 
inflation ||l|,^ , connected to the existence of fractal basin 
boundaries in the set V. For a given radius R of 2?, it 
is always possible to find a gap of energy 6 E* (which 
depends on R^) such that for Eq e 5 E* , the exit to in- 
flation is chaotic, namely, small fluctuations (of order of, 
or smaller than R) of the initial conditions will change the 
asymptotic state of the orbits from collapse into escape 
to inflation, and vice versa, after visiting the neighbor- 
hood of the saddle-center. The chaotic exit to inflation 
may be characterized in terms of the partition of energy 
(2): for a particular orbit, the partition into the approxi- 
mately conserved energies E^yp and Erot about the criti- 
cal point will determine the asymptotic state of the orbit, 
namely, collapse or escape to inflation, whether respec- 
tively Ehyp < or Ehyp > 0. However, we are no longer 
able to foretell precisely what will be the value of Ehyp 
endowed to a given orbit of 2? when it visits a neighbor- 
hood of E, due to the high sensitivity of Ehyp with the 
fluctuations of the initial conditions, a consequence of 
the chaotic nature of the set P. In what follows, we will 
use standard statistical methods to examine this instabil- 
ity of the partition of the energy, a procedure which will 
lead us into the recognition of new statistical patterns for 
chaotic sets in the phase space. 

To start, let us make a plot of Ehyp versus the ordered 
values of one of the canonical variables, say a, especi- 
fying the initial conditions in the set P, for R m 10^*^, 
m — 2 and {Ecrit ^ Eq) = 10~^°. We fixed the value of 
A = 3/2, such that uq = 1 and Ecrit = 3/2. Fig. 1 cor- 
responds to 3, 031 initial conditions, sampled randomly 
and uniformly in P, and the points were connected for 
successive a's to a better visualization. The result shows 
a highly disorganized " signal" , displaying structures in 
all scales of a. Analogous highly disorganized "signals" 



appears for other variables (/?, Pa, or p^. Also, if another 
sample of points is taken randomly and uniformly in P, 
we obtain a "signal" with the same general aspect but 
completely different in all the details, appearing not to 
have any connection with the previous figure. However 
a pattern of the signal is reproducible, which is its his- 
togram. 

Our basic procedure is to realize the following numer- 
ical experiments with the system: from a given ensem- 
ble V of initial conditions, we generate orbits of the sys- 
tem by the full Hamiltonian dynamical equations derived 
from (1). These orbits visit a linear neighborhood of the 
saddle-center E. In this linear neighborhood, we cal- 
culate the approximately conserved quantity Ehyp asso- 
ciated with each orbit ; for practical purposes, Ehyp is 
evaluated at the point of closest approach to E, in the 
sense of the Euclidean measure. We construct the his- 
togram function Ps^y^ = A N (Ehyp)/ A Ehyp as usual: 
the Ehyp axis is divided into a large number of equal small 
bins A Ehyp and the histogram is given by the number 
of orbits assigned with hyperbolic energy in the interval 
[£;gp- A Ehyp/2, e\11^+A Ehypl2\ , when the orbit visits 

the neighborhood of E. Here -E^y^ is the value of the hy- 
perbolic energy centered in the i-th bin. We fix the origin 
of the histogram in E^^^ , which is the minimum value of 
Ehyp obtained in the numerical experiment; according to 
(2), E;™" > -{Ecr^t - Eo). We also note that due to (2) 
the histogram of the signal associated with Erot results 
merely from a shift of the histogram for Ehyp. 

A careful examination shows us that the form of the 
histogram for smaller values of Ehyp is sensitive to the 
geometry of the sampling set V. We must therefore es- 
tablish a criterion to fix the geometry of V such that no 
bias due to the geometry of the sampling volume man- 
ifests itself in the histogram, and such that the statis- 
tical distribution function obtained from the histogram 
will be connected only to the non-intcgrability of the sys- 
tem. We establish that the domain V of randomly and 
uniformly distributed phase space points must have a ge- 
ometry such that it yields histograms A N {Ehyp)/ A Ehyp 
or A N {Erot) / A Erot for the integrable Hamiltonian H — 
i (p^ + rn^ f '^) + j2 Pa ~ ^ o.'^ having the form of a plateau 
(cf. also the caption of Figs. 2). In the case of the non- 
integrable Hamiltonian (1) the histogram constructed 
from the same set of initial conditions V deviates from 
a plateau, reflecting the nonintegrability of the system. 
Figs. 2(a,b) show 30,000 initial conditions sampled ran- 
domly and uniformly in V constructed for m = 10, 
R w 10~^ and {Ecnt — Eq) = 10^^", projected on the 
planes (95, Pi^) and {a,pa). Fig. 3 shows the histogram 
generated from the initial conditions set T> of Fig. 2 by 
the non- integrable Hamiltonian (1). 

At this point our main question is: what are the possi- 
ble candidates for the statistical energy distribution law 
associated with the histogram function? The present his- 
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tograms have an excess in the tails, namely, in the region 
of large positive hyperbolic energies, relative to the cor- 
responding Maxwell-Boltzmann distribution. Indeed, a 
trivial log-linear plot of the histogram shows a non linear 
form for the tails, exhibiting an anomalous behaviour in 
contrast to the case of the normal Maxwell-Boltzmann 
distribution. This is not surprizing since we are dealing 
with a system were long range interactions are present, 
and where relevant phase-space sets are chaotic and have 
a fractal-like structure. The system is certainly not in the 
realm of Boltzmann-Gibbs statistical mechanics, which 
essentially focus on non-interacting or short-range inter- 
acting many-body Hamiltonian systems. After extensive 
numerical analysis, we obtained that all histograms ex- 
amined are fitted by the energy distribution function 

P{x)^CoP{f3xy\ (3) 

where Cq is an adimensional normalization constant and 
X = {Ehyp — E™yp). The parameters of the distribution 
have values 7 = 0.6354 ± 0.02527, CqP^"'^ x 10^°-^ = 
12.496 ± 0.5108. These are the values for the best fit 
of the histogram of Fig. 3, with the square of the vari- 
ance = 0.06022, obtained using the PAW packet from 
CERN. Within statistical errors, the value 7 = 0.6354 
corresponds to the best fit of all histograms generated 
with the dynamics of (1). For comparison purposes, all 
histograms were normalized by fixing the peak with the 
value 100. 

The distribution law (3) produces an accurate fit for 
almost the whole range of energies of the histograms. 
The continuous line in Fig. 3 shows the best fit of (3) to 
the histogram. The whole histogram is composed of 100 
bins, each with A Ehyp = 7.33993555 x lO^^^. To obtain 
the fit with the value — 0.06022, we have neglected 
the extreme of the tail of the histogram, more specifically 
from the 55*'' bin on. If we have extended the solid curve 
beyond this bin, a discrepancy between the curve and the 
histogram would appear. This discrepancy at the end of 
the tail has the following explanation. The cut-off at the 
largest values of E^yp is due to the finite volume of the 
sampling set T>. If we increase R (consequently increas- 
ing the volume) and maintain the same uniform density 
of points, the cut-off is shifted to the right and the tail 
of the histogram raises improving the fit in this region. 
The left extremity is a peak. If we extended (3) to values 
of X very close to zero, the peak might suggest a diver- 
gence in the density of states - a property which is not 
physically reasonable for statistical ensembles - but this 
is not the case. Indeed the smallest values of E^yp in 
our histograms are associated to points that are located 
at the boundaries in (a,pa) of the linear neighborhood 
of the saddle-center E, and the evidence from the his- 
tograms is that, as we increase the volume of V, there is 
a saturation in the corresponding increase of the peak of 
the histograms. Therefore the distribution function (3) 



must surely be corrected, to be finite in the neighbor- 
hood of X = 0. The saturation is possibly due to the fact 
that, as we enlarge the volume, we approach the limit 
of validity of the linear approximation, and higher order 
terms in the Hamiltonian may be contributing already to 
make the density of states finite there. The parameter 7 
is adimensional while 1/(3 must have the dimension of en- 
ergy. We may tentatively interpret the latter as defining 
a temperature for the saddle-center. For the integrable 
case, 7 = 0. 

The following characterist properties of the histograms 
were verified numerically, through exhaustive numerical 
experiments realized with randomly and uniformly dis- 
tributed distinct sampling sets P having the same geom- 
etry described in Figs. 2: 

(j) the form of the histogram is scale invariant, that is, 
it is independent of the characteristic radius R of the 
sampling set T) of initial conditions, modulo a change of 
bin; 

(m) the form of the histogram is independent of the pa- 
rameter {Ecrit ^ Eq) that allows to characterize if the 
orbits are dominantly collapsing or escaping, after visit- 
ing the neighborhood of the saddle-center E; 
[Hi) the form of the histogram is independent of the mass 
parameter m; 

{iv) the histogram has the same form also for a model 
with dust, what lead us to conjecture its independence of 
the particular equation of state p — \ p, with < A < 1; 
[v) within statistical errors, the value 7 = 0.6354 corre- 
sponds to the best fit of all histograms associated with 
the Hamiltonian (1). 

We were not able to identify (3) with any known sta- 
tistical distribution. However, we are aware that this be- 
haviour is typical of tails of distributions arising in Tsal- 
lis generalized statistics that were proposed in order 
to accomodate, at least in part, systems which have an 
anomalous behaviour with respect to Boltzmann-Gibbs 
statistics. 

We now exhibit numerical evidence that the constant 
7 appearing in (3) may characterize the fractal nature 
of the sets 2?, more specifically, it may give a measure of 
the fractal dimension of basin boundaries in these sets. In 
fact, by using a box counting method ||l0|] with the un- 
certainty code defined by collapse/escape, we obtained 
the uncertainty exponent a = 0.36931 for the set V with 
TO = 10, iEc„t-Eo) = 10-1" and R = 10"^ (cf. Fig. 5). 
Within statistical errors, this same value was obtained 
for all initial conditions sets with Eq E S E*, used to con- 
struct histograms in our numerical experiments, implying 
that 

7 + awl. (4) 

As a consequence 7 appears to be directly related to the 
fractal dimension d of the basin boundaries of 2?, since 
in our boxcounting calculation d = N — a, where = 3 
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is the dimensionality of T>. The above restriction to sets 
V with Eq & 6 E* is due to the fact that in these sets we 
may define the uncertainty code collapse/escape; for sets 
with Eq ^ S E* , other uncertainty code may be used. 
We have also examined the Hamiltonian 
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2AB^-Eo = 0, 



(5) 



describing the dynamics of anisotropic spatially homo- 
geneous Bianchi IX cosmological models, with scale fac- 
tors A and B, plus a cosmological constant and a perfect 
fluid in the form of dust. This system is chaotic and 
has a saddle-center in the finite region of phase space [D . 
We made the same construction of histograms for hyper- 
bolic/rotational energies in a linear neighborhood of the 
saddle-center; in this case, we have the burden of ob- 
taining Moser's normal coordinates in order to calculate 
properly the hyperbolic/rotational energies associated to 
each orbit when visiting the neighborhood of the saddle- 
center. Fig. 4 shows the histogram generated from 30,000 
initial conditions. The geometry of the initial conditions 
set satisfies the same criteria of the previous case (1). The 
continuous line corresponds to the best fit of (3), with 
7 = 0.71458±0.04060, Co P^''^ x IO^Ot = 11.601±0.7711, 
and = 0.06685. The histograms consists of 65 bins, 
each with value AEhyp = 3.58919436 x 10"^^ The best 
fit was realized in the interval of the 2"'^ to the 30*'' 
bin. Within statistical errors, all histograms associated 
with (5) are fitted by the distribution function (3) with 
7 = 0.714. Properties (i) and (m) were also verified for 
the present case. A boxcounting for the initial condi- 
tions set yields the uncertainty exponent a = 0.28103, 
confirming the relation 7 -I- a = 1. 

We end this paper by making the conjecture that the 
statistical distribution (3) is an universal characteristic of 
saddle-center critical points, independent of the particu- 
lar dynamical system in consideration. We are presently 
working in this direction by making the same constructo 
for two other Hamiltonian systems presenting a saddle- 
center in their phase space, not in the realm of inflation- 
ary cosmologies: the Henon-Heiles system jllj and the 
chemical system examined in Ref. |T3]. Our preliminary 
results indicate the same distribution function (3) with 
distinct values of 7 for each system. We also conjecture 
whether other critical points of different nature can be 
universally characterized in this way. 

The above results may have some interesting physi- 
cal applications for the dynamics of preinflationary mod- 
els immediately before the universe enters the inflation- 
ary regime or recollapses, both alternatives depending 
on the typical fluctuations in the system. For instance, 
we may use the statistical function (3) to calculate the 
average value of the hyperbolic energy for ensembles 
of universes. As we know [|l],|j the sign of this av- 
erage hyperbolic energy determines the "average" long 



time behaviour of the ensemble, that is, collapse or es- 
cape to inflation. A straightforward calculation shows 
that the sign of the average energy depends basically on 
((1 - 7) /? E'^yl'' + (3 EJ^'y'^), where EJ^yl'^ is the maximum 
value of the hyperbolic energy in the histogram. The rel- 
ative value of E™yp decreases linearly (and may exceed 
in absolute value E^y^) with the increase of {Ecrit — Eq), 
implying that larger/smaller values of {Ecrit — Eq) deter- 
mine the dominance of collapse/escape in the dynamics 
of the ensemble T> of universes. 
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FIG. 1. Signal resulting from the distribution of hyperbolic 
energies Ehyp versus the ordered initial values of a. We consid- 
ered 1,035 initial conditions sampled randomly and uniformly. 
Points were connected for successive a's, for better visualiza- 
tion. 
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FIG. 2. Set V of 30,000 initial conditions, uniformly and 
randomly distributed, for the case m = 10, Ecrit — Eq — 10^^" 
and R = 10~*, projected in the planes (fi,p<p (Fig. 2a) and 
{a,pa){Fig.2b). The particular geometry taken for the initial 
conditions set is such that (i) it satisfies the Hamiltonian con- 
straint (1) and, (ii) the associated histograms for the integrable 
case have a form of a plateau; this demands that Fig. 2a be 
an ellipse with axes having the ratio m. 
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FIG. 3. Histogram of the TMstMbution of Ehyp, associated 
with the 30,000 initial conditions of Fig. 2. The continuous line 
is the best fit of the distribution (3) with 7 = 0.6354± 0.02527, 
Co /3^"^ X 10^°^ = 12.496 ± 0.5108 and = 0.06022. The 
horizontal axis is scaled with the factor 10~^°. 
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FIG. 4. Histogram of the distribution of Ehyp, generated 
from 30,000 initial conditions for the anisotropic case, with 
A = 1/4, Ecru - Eo = 0.5 X 10"^" and R = lO"''. The 
continuous line is the best fit of (3) with 7 — 0.7146 ± 0.0406, 
Co (3^-'' X 10^°T = 11.601 ± 0.7711 and = 0.06685. The 
horizontal axis is scaled with the factor 10~^°. 
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FIG. 5. Plot of / = feo in a log-log scale, where e is the un- 
certainty radius about 20,000 initial conditions taken inside the 
set D, and / is the fraction of T> of uncertain initial conditions 
with the uncertainty code collapse/escape. We considered 20 
values of e inside the interval [0.9 x 10"*, 1.8 x 10"**]. The best 
fit renders a « 0.36931. 
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